#ifndef PLANT_SPDE2D_H
#define PLANT_SPDE2D_H

#include "../Main/Signal.h"
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>

const double SQRT_TWO = 1.414213562373095;

class plant_SPDE2D:public block
{
	int    mode;
	int    LatticeSize;
	double *Km2[4];

	// States
	double *z1;
	double *z2;
	double *z3;
	double *z4;
	
	double T;
	double W;
	double c;
	double c2;
	double sigma;
	double **lambda;

	double dt_con; // dt_con is the dt used in numerical integration. 
	               // block::dt is the time interval for output.
	double sqrt_dt_con;
	double dt_conservative_factor;

	// Random number generator supported by GSL
	const gsl_rng_type * noise_type;
	gsl_rng * noise;

	double eq_c(double T,double W);
	double eq_c2(double T,double W);
	double eq_sigma2(double T,double W);

	inline int index2D(int m,int n){return m*(mode+1)+n;};
	
	inline double phi1mn(int m,int n,int ix,int iy){
		return 2.0*sin(2*m*pi/LatticeSize*ix)*sin(2*n*pi/LatticeSize*iy)/pi;
	};
	//double phi2mn(int m,int n,int ix,int iy);
	inline double phi2mn(int m,int n,int ix,int iy){
		double coeff;
		if((m == 0 || m*2==LatticeSize)&&(n==0 || n*2==LatticeSize)){
			coeff = 1.0/pi;
		}
		else if(m!=0 && 2*m!=LatticeSize && n!=0 && 2*n!=LatticeSize){
			coeff = 2.0/pi;
		}
		else{
			coeff = SQRT_TWO/pi;
		}
		return coeff*cos(2*m*pi/LatticeSize*ix)*cos(2*n*pi/LatticeSize*iy);
	};

	inline double phi3mn(int m,int n,int ix,int iy){
		double coeff;
		if(n==0||n*2 == LatticeSize){
			coeff = SQRT_TWO/pi;
		}else{
			coeff = 2.0/pi;
		}
		return coeff*sin(2*m*pi/LatticeSize*ix)*cos(2*n*pi/LatticeSize*iy);
	};

	inline double phi4mn(int m,int n,int ix,int iy){
		double coeff;
		if(m==0||2*m==LatticeSize){
			coeff = SQRT_TWO/pi;
		}else{
			coeff = 2.0/pi;
		}
		return coeff*cos(2*m*pi/LatticeSize*ix)*sin(2*n*pi/LatticeSize*iy);
	};
public:
	plant_SPDE2D(int M,int L,double idt);
	~plant_SPDE2D();
	virtual void update(double sysTime);
	virtual void reset();
	void update(double sysTime,double T,double W);
	void output(string filename);
	double getM2();
	double getR2();
};

void h2mode_2D(double *h,int L,double *z,int M);
void mode2h_2D(double *z,int M,double *h,int L);
double randn();
#endif